area1=c(15,14,21,13,9,10)
area2=c(4,6,10,9,8,5)
area3=c(10,11,14,10,9,15)
area4=c(8,3,4,1,2,5)
trees=data.frame(area1,area2, area3, area4)
trees

Mi=c(12,15,14,21)
mi=c(6,6,6,6) 
N=10 
M=sum(Mi) 
m=sum(mi) 
n=4 
Mbar=M/N
area.means=colMeans(trees, na.rm=TRUE)
area.means
muhat=(1/Mbar)*sum(Mi*area.means)/n 
s2b=sum((Mi*area.means-Mbar*muhat)^2)/(n-1)
s2i=apply(trees,2, var, na.rm=T) 
vhatmu=(1-n/N)*(1/(n*Mbar^2))*s2b+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(s2i/mi)) bound=2*sqrt(vhatmu)

tauhat=muhat*M 
vhattau=vhatmu*M^2 
boundtau=2*sqrt(vhattau)

rbind(Mi,mi,N,M,m,n,Mbar,muhat,s2b,vhatmu,bound,tauhat,vhattau,boundtau)